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Abstract 

We study a class of 2-dimensional Hamiltonian systems H{x,y,px,Py) = \{p'i 
+ V{x^y) in which the plane x=Px=^ is invariant under the Hamiltonian flow, so 
that straight-line librations along the y axis exist, and we also consider perturbations 
5H = 5 ■ F{x,y,px,Py) that preserve these librations. We describe a procedure 
for the analytical calculation of partial derivatives of the Poincare map. These 
partial derivatives can be used to predict the bifurcation behavior of the libration, 
in particular to distinguish between transcritical and fork-like bifurcations, as was 
mathematically investigated in [1] and numerically studied in [2]. 

1 Introduction 

We study 2-dimensional Hamiltonian systems 

H{x,y,Px,Py) = lipl + pl) + V{x,y) (1) 
with a potential satisfying 

dV 

— (0,y) = for all yeR. (2) 

Then the {y , Py)-pla.ne x = Px = is invariant under the Hamiltonian flow and 
thus the system will librate on the y-axis. Choose one of the libration families, 
parametrized hy e := E — Eq & I for some fixed energy Eq and a suitable open 
interval /. The family consists of closed orbits jeit) = (0, t), 0,pj;(£, t)), with 
Py{e,t) = y{e,t)), which we let start, say, at their maximal value of y, that is at the 
point (0, ?/inax(^)7 0, 0) G M^. Let T{e) > denote the period of 7^. 

For each e we may use = as a Poincare surface of section PSS at the starting 
point of the orbit, and the (a;,Pa;)-plane as the PPSS, the projected Poincare surface 
of section. As the canonical coordinates q and p in the PPSS, we may choose x and 
Px- Then the Poincare map defines a symplectic family 

Q = Q{(i,p,£) /ON 

P = P{q,p,e) 

on an open neighborhood of the e-interval x x J in the (g,p, 5)-space M^. Note 
that 

A:=OxOxI = {{0,0,e)\eeI} (4) 

itself is a fixed point branch of this family, and we propose to study the bifurcations 
that may occur along this branch. 
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Among the questions we ask about these bifurcations, there is one concerning the 
behavior of a bifurcation under a small deformation of the Hamiltonian. Let S denote 
a small deformation parameter and let us add a deformation term S ■ F{x,y,p^,py) 
to the Hamiltonian: 

H{x,y,p^,Py,6) = \{pI+pI) + V{x,y) + 5F{x,y,p^,Py), (5) 

and let F satisfy the following 'libration preserving condition': 

dF OF 

— (0,2/,0,p,) =0 and _(0, 0,p,) = 0. (6) 

Then for fixed 6, the {y,py)-p\a.ne x = px = will still be invariant under the 
Hamiltonian flow, so the system will still librate on the y-axis, and we ask how the 
bifurcation behavior found at 5 = will change if we turn on the parameter 6. 

The symplectic family defined by the Poincare map now depends on two param- 
eters e and S: 

Q = Q{q,p,£,S) 

P = P{q,p,e,6). ^'^ 

Think of an arbitrary Sq being chosen. We will ask if the fixed point (0,0), of the 
undeformed system, is singular at e = Sq and if so, what are the properties of the 
bifurcation and their behavior under deformations 6^0. Using the answers to 
these questions could be read from those 38 partial derivatives up to third order of 
P and Q that involve the parameters e and 6 at most in first order, at the single 
point (OjOjSoiO) — if we only knew them. The purpose of the present note is to 
describe a procedure for the calculation of these partial derivatives of the Poincare 
map from the potential V{x,y) and the deformation term F{x,y,Px,Py). 



2 Numerical prerequisites 

To start the procedure, for a given Sq, we will need to know first of all the closed 
orbit itself, 7eo(t) = (0, y(t, Eq), 0,Py{t, Eq)), that is we have to solve the equation 

dV 

2/+^(0,^/) = (8) 

to the initial condition y{0) = ymaxi^o) and y{0) = 0. The value 2/max(^o) satisfies 

V{0,ymU^o)) = Eo + Eo. (9) 

Within the chosen domain of libration it will be the larger of the two solutions of 
this equation and can be determined that way. The function y{t) := y{t, Eq) will be 
periodic of a period T{eo) > 0. 

Furthermore, we will have to solve the linearized Hamiltonian equation along 
this closed orbit, that is we have to know the fundamental system (^1,^2) of the 
linear equation 

e + ^(o,2/(t))e = o (10) 
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to the initial condition 

/ fifO) £,(0) \ / 1 \ 
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6(0) 6(0) 
6(0) 6(0) 











^ 





as well as the fundamental system (?]i,772) of 



fi+^(0,y(t))rj^O (12) 



to the initial condition 



^i(O) r/2(0) 
?7i(0) 7)2(0) 




(13) 



To say that these five functions yit), £,i{t), 6(^)) Vii't)^ Viif) and their first derivatives 
must be 'known' means that they are known numerically on the whole period interval 
[0,T(£o)]- A computer program implementing the procedure for the calculation of 
the partial derivatives of the Poincare map at (0, 0, £q, 0) will have to treat them as 
known functions. But beyond that no further differential equations will have to be 
solved. 



3 The Poincare map 

The Poincare map is produced by the Hamiltonian flow. To facilitate the handling of 

higher partial derivatives, wc will gradually shift from natural 'speaking' notations 
like ( x,y,Px,Py) to a simple enumeration of variables by upper indices. Partial 
derivatives by these variables will then be denoted by corresponding lower indices. 
We begin by writing 



a 

4 



1 -^x 



a 

p4 



= y 



for the independent variables in M and 

= x*(^,a^a^a^a^^) (15) 

with i — 1, ... ,4, for the components of the Hamiltonian flow at a fixed 5, with 
initial conditions a^, . . . , : 

x%0,a\a^,a^,a^,5) = a\ (16) 

For fixed e and 6, the starting point in the PSS corresponding to a given point {q,p) 
in the PPSS is described by 

— q 

= p ^ ^ 
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where the y-component is defined imphcitly by 

+ Viq, yiq,p, e,5)) + 5F{q, y{q,p, e,5),p,0)^e + Eo (18) 

and j/(0,0,£, 0) = ymax(£)- From this starting point, the flow fine will travel for a 
time T — T{q,p, e,S) > until it hits the PSS Py — again, so implicitly this time 
is given by 

x^iT{q,P,s,S),q,y{q,p,e,S),p,0,S) = 0, (19) 

in the notation (15) of the flow, and by the reference condition ^(0, 0, £, 0) = T{£), 
the period of %. The Poincare map can now be described as 

Q{(1,P,£,S) = x\T{q,p,e,S),q,y{q,p,e,5),p,0,5) .^g) 
P{q,P^£^S) = x^{T{q,p,e,5),q,y{q,p,£,5),p,0,S). 

Taking partial derivatives, we obtain the the partial derivatives of the Poincare map 
in terms of partial derivatives of the Hamiltonian flow and of partial derivatives of 
the flow time function T{q,p, e, 6) and the starting point function y{q,p, s, 6). That's 
what we do next. 



4 Taking derivatives of the Poincare map 



We are now unifying the notation of the independent variables of the flow, currently 
written as {t,a^,a'^,a^,a'^,6), to (a°, . . . , a^). Derivatives are denoted by lower in- 
dices, so for instance would mean 



^02 — 



dx^ 



dtd{a^) d{a^ 



(21) 



and so on. Next, we write the two components Q and P of the Poincare map as 
compositions 

Q^x^oZ 
P^x^oZ, 



(22) 



or Q{q,p, e, 5) = x^{Z^{q,p, e,5), . . . , Z^{q,p, e, 5)) and analogously for P, where the 
six components Z^, . . . ,Z^ are given, according to (20), by 



(23) 



Finally, let us also enumerate the independent coordinates {q,p, e, S) in the product 
of the PPSS with the parameter plane by writing 



Z\q,p,e,S) 


T(q,p,s,S) 


Z\q,p,e,5) 
Z^{q,p,£,6) 
Z^{q,p,e,6) 
Z^{q,p,e,5) 


:= q 


■= y{q,p,£,^) 


■.= p 


:= 


Z^q,p,e,5) 


:= 5. 



u 



Q 
P 

e 



(24) 
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again denoting partial derivatives by lower indices, like in Z| = = I shall 
try to be consistent in using greek letters A, . . . for the a-related indices that run 
from to 5 and use roman letters m, ... = 1,2,3,4 for indices refering to the 
M- variables. 

The partial derivatives we are after, like P^^g = (.t'^ o Z)ii3, can now neatly be 
written in terms of partial derivatives of flow, flow time and starting point functions 
as 

{x'oZ)e = x\Z^ 

ix^oZ)em = x\^Z^Zi:^ + x\Z^^ (25) 

{x' O Z)emn = ^XiiuZ^ZI^Z^ + x\^{Z^Z!f^^ + ^4^n€ + ^n^tm) + ^\^£mn- 

That's just the chain rule, so it holds everywhere. But for our application, we 
need to know the left hand sides only at the point (w\ u^, u^, u^) = (0, 0, Eq, 0), and 
so on the right hand side we want to know the Z^ at (0, 0, Eq, 0) and the x% at 
{T{Eo),0,y^^{Eo), 0,0,0). 



5 Derivatives of the starting point function 

First we will take care of the derivatives Z^. For A = 1, 3, 4, 5 this is easy, because 
for those A the Z^ are given by 

Z} = Zl = Z\ = \ and 

(26) 

for all other lower indices *, in particular for all derivatives of order > 2. It remains 
to determine the derivatives Z\ and Z\ of Z^ = T{u^, . . . u"^) and Z"^ — y{n}, . . . u^). 
In the present section we will calculate Z| up to second order. 

The defining equation (18) for the starting point function y{v}, . . . ,u^) in the 
li-notation becomes 

■ + V{u\ y{u\ u")) + u^F{u\ y{u\ u'),u', 0) = + Eq. (27) 

As a first step, we will determine the eight derivatives 

yi,2/2,y3 and 2/11,2/12,2/22,^13,2/23 (28) 

at (0,0, £0,0). Since here the variable — b \% not involved, we may put — 0. 
Taking first derivatives by m)-, , gives 

Vi + V2yi = 

+ 1^22/2 = (29) 
^22/3 = 1 
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at all {u^ , , u"^ , 0) . Differentiating further, we get 

Vu + Vuyi + {V21 + V22yi)yi + V2yn = 
Vi2yi + V22y2yi + ^^2^/12 = 

1 + V^22y2Z/2 + ^"2^/22 = (30) 

{V21 + V22?/l)Z/3 + V2y31 = 

V22?/22/3 + V2y32 = 0, 

again at all {u^, u^, u^, 0). Specializing now to (0, 0, Eq, 0) and using our assumption 
|j^(0, y) = on the potential, which implies 

V,{0,y) = V,2{0,y) = (31) 

for all y, we obtain 

yi = ?/2 = yi2 = Z/13 = ?/23 = 

and 1/3 = 



2/11 



1 

V2 

Vn (32) 

V2 

1 

y22 = -77 
1/2 

at (0, 0, £0) 0) as some first 'end results' on the derivatives of the starting point func- 
tion in terms of the derivatives V2{0,y^g^{so)) and Vii(0, ?/max(^o)) of the potential. 

Let us now consider the derivatives ^4, yu, and y24 involving the deformation 
parameter u'^ — 5. Differentiating (27) first by and then in addition by resp. 

we obtain 

V2y4 + F = 
(V^2i + 1422/1)2/4 + ^2/41 + i^i + i^22/i = (33) 

^222/22/4 + ^2/42 + ^22/2 = 

at all {u^ , , , u"^) . Specializing again to (0,0,60,0) and using (31), (32) and the 
assumption (6) about F{x,y,Px,Py), we get 

2/14 = 2/24 = 

i^(0,2/max(eo),0,0) (34) 



and 1/4 



V'2(0,ymax) 



Note that with (26), (32) and (34) wc have determined all partial derivatives 
and at (0,0,£0)0) for A > 1 and at most one of the indices £ and m being 3 
or 4. For convenience, let us collect these results. First for the 'easy' A's. Here the 
first derivatives are, everywhere: 





e=i 


£=2 


£=3 


£=4 


A=l 


1 











A=3 





1 








A=4 














A=5 











1 



(35) 
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The higher derivatives therefore are zero, in particular = for A = 1,3,4,5. 
For A = 2 we have found at {u^, v? , u^, u'^) = (0, 0, Eq, 0): 

Zf = 

zl = 



^3 



1 

F 

2 _ ^11 



72 
^22 



2 



^?2 = 



1 



zh = 

ZI3 = 

= 

Z24 = 0' 



(36) 



where the derivatives of V have to be taken at {x,y) = (0, ?/max(£o))) the value of 
F at {x, y,Px,Py) = (0, ymax(£o), 0, 0). In the next section, we turn to the remaining 
case A = 0. 



6 Derivatives of the flow time function 

Although we are interested in the partial derivatives of x* o Z for 2 = 1,3 only, we 
will also become involved with the Pj^-component a;*^(a°, . . . , a^) of the flow, because 
the defining condition (19) of the flow time function Z° is 

x^oZ = 0. (37) 

This is also the reason, by the way, why the linear equation (12) will come up in the 
calculations. — From (37), using the first equation of (25) for i = 4, we obtain 

^oZe — —^iZ] — x^Zl — x\Z\ — x\Z\ — x\Z\. (38) 

The Z^ on the right hand side arc known at (0,0,£0)0). Correspondingly, the x\^ 
for A = 0, . . . , 5 are meant to be taken at (a°, . . . , a^) = (T(£o), 0, ymax(^o)) 0) 0) 0)- 
We still have to determine them, but we certainly know x\ there, since by the 
Hamiltonian equations 

dV 

xtit, 0, yu^U^o), 0, 0, 0) = py{t) = -—(0, y{t, £0)) (39) 

along the closed orbit 7£,j, and therefore 

x^(T(£o),0,y„ax(£o), 0,0,0) = -^2(0, yn.a.(£o))- (40) 
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So from (35), (36) and (38) we get 



70 _ 1 ^4 



Z2 



\ (41) 

-1- A 



7O _ ™4 I ^ 4 



Similarly, from the second equation of (25) for i = 4 we now have 

4 2 4 A . ^''^ 

at (0, 0, £0) 0), with the to be taken at {T{eo), 0, 2/max(£o); 0, 0, 0), as before. And 
here we leave it for now: the Z^ on the right hand side, at (0, 0, £oi 0), are all known 
from (35), (36) and (41). A computer will understand (42) as given, and for us, 
there is no point in writing out the formula in great length before we know more 
about the x'l^, in particular before we know which of them will vanish anyway. 

The same reasoning applies to the last Z^ that are still missing, namely the third 
derivatives and Z^^^ of the starting point and flow time functions, which might 
be needed in the computation of the third derivatives of the Poincare map. In fact 
they will not be needed, because they enter the third equation of (25) for i = 1,3 
with coefficients x\ and Xq, which will soon be seen to vanish for i = 1,3. This is 
just one of the details of the problem to calculate all the x% (T(£o) , 0, j/max(£o) ) 0, 0, 0) 
that we need. To this problem our bifurcation analysis is now reduced and it will 
be solved in the remaining sections. 



7 Equations for the flow derivatives 

The four components of the Hamiltonian flow, see (15) and (16), are written cur- 
rently as 

x' ^x'{a°,...,a^) for i = l,2,3,4 (43) 

with a° denoting the time t and — S, while (a^, . . . , a^) is the initial point. As a 
bookkeeping device we now introduce a 'fifth fiow component' by 

x'(a°,...,a^) :=a^ (44) 

But we also use x^, . . . ,x^ as the names of the independent variables of the Hamil- 



9 



tonian vector field v = (f \ . . . , w"^), which is then given by 

v^{x^, . . .,x^) 
v^{x^, . . . ,x^) 
v'^{x^, . . . , x^) 

The Hamilton equations become i;^(a) = v'^{x{a)) for r=l,2,3,4, and as in (25) we 
obtain 

xl = v^x\ 

= <j4^^M + <^AM (46) 

/y>T* 'I /y^ fyJ I T'J I /y*?' i-y»«? I rr*^ rin3 \ I n\'f /yt 

The summation indices A; run from 1 to 5, while any Xj/J-jV G {0,1,2,3,4,5} 
are admitted. The equations hold everywhere, that is the flow components and 
their derivatives may be taken at any a := (a°, . . . , a^), the vector field components 
and their derivatives then at the corresponding x(a)) = {x^{d), . . . ,x^{d)). Also 
note that the time derivative of the flow is the partial derivative by a°, so on the 
left hand sides we might have written x^^, x^^^, Xq^^^, instead of x\, x\^, Xx^v- 

To apply (25), we will only need to know the x%{T{eo), 0, ?/max(£o), 0, 0, 0). But 
in order to determine these numbers, we will also have to consider the functions 
a;!|,(t, 0, |/max(£o)j 0, 0, 0) on the interval [0,T(£o)], for which we now introduce the 
notation 

4(t) := 4(i,0,y^ax(£o), 0,0,0). (47) 
Correspondingly, we write 

vl{t) := vl{x{t,0,y^^{so), 0,0,0)). (48) 

Then from (47) we get 

xlit) = vlit)x\{t) 

= vl^{tK{t)x^^{t) + vl{t)x\^{t) 
i^x,At) = vl^,{t)x\{t)xlit)xtit) (49) 

+vl{t)x\,M) 

from which we will now proceed to determine the functions (47). We only consider 
r = 1, . . . , 4 as there is no need to write equations for x%{t), since of course 

xl = l and xl{t) = (50) 

for all other partial derivatives of x^{d) — a^. 



x^ + x^Fs{x\...,x^) 
x^ + x^F4(x\...,x'^) 
-Vi{x\x'^)-x^Fi(x\. 

-V2{x^,x'^) - x^F2{x^, . 



.,x') 
.,x'). 



(45) 
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8 Calculation of the first order flow derivatives 



Note that all the {t) on the right hand sides of (49) are known functions in the 
sense agreed upon in section 2 on the numerical prerequisites, since 

x{t, 0, yr.U^o),0, 0, 0) = (0, y{t),0, y{t), 0), (51) 

which is to be used in (48). In particular, let us tabulate the first derivatives vl{t): 





i=l 


i=2 


1=3 


i=4 


i=5 


r=l 








1 








r=2 











1 


m 


r=3 


-Vuit) 














r=4 





-V22{t) 








-F,{t) 



(52) 



Here we write F4{t) :— -F4(0, y{t),0, y{t)) and Vu{t) :— Vii(0, y{t)) and so on, in line 
with the notation introduced in (47) and (48). Note that Fi{t) = F^lt) = by our 
assumption (6) on F{x,y,px,Py). As a first consequence of (49), we see that the 
x\{t) for r, A e {1,2,3,4} satisfy the homogeneous linear differential equations 



XX 



+ Vn{t)xl 



and 



X\ 



+ 











(53) 



(54) 



V22{t)xl 

These are just the first order systems corresponding to (10) and (12), and the initial 
conditions, as we see from (16), are 



x\{0) = 1 if r = A, 
a;^(0) = if r ^ A. 



and 



(55) 



But this shows that the x\ (t) for r, A 
numerical prerequisites, more precisely 

x\it) xlit) 
xl{t) xlit) 



,4 are known functions from the 



and 



while 



xlit) xlit) 

xm xiit) 



Ut) Ut) 

viit) mit) 
viit) mit) 



(56) 



(57) 



(58) 



What about A 



xlit) = for r + A odd, AG {1,2,3,4}. 
and A = 5? We know x^it) = x^it) = v^it), so from (45) we have 

xlit) = 

xlit) = yit) 

xlit) = 

xlit) = -V^it). 



(59) 
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Looking at A = 5, we see from (47) and (52) that (53) is satisfied also in this case, 
but (54) has to be replaced by the inhomogeneous system 



Xj + V22{t)x 



-F2{t). 



(60) 



The initial conditions are ^5(0) = 0, as we see again from (16). In particular we 
have 

xlit) = xlit) = 0, (61) 

and since we have the fundamental matrix (57) of the homogeneous system (54), we 
obtain the the remaining two functions xHt) and x^t) by variation of constants as 




xiit) xtit) 




4i 



r) 
-4{r) 




• (62) 



Note that with (50) and (56)-(62) all flow derivatives x\it) of first order are now 
determined. Among them the x^t) and x^it) have been seen to vanish for i = 1,3, 
as announced at the end of section 6, and so there is in fact no need to determine 
third derivatives Z^^^^ and Z^^^ of the starting point and flow time functions. 



9 Calculation of higher order flow derivatives 



Just as the first equation of (49) led to (53) and (54), so the other equations of (49) 
show that the higher order flow derivatives x%it) satisfy differential equations 



X<i 



xl + Vnit)xl 



and 



with 



X<i 



X-i. — 



xi + V^22(t)a;| 



gl{t) 

gl{t) 
g%it), 



gUt) 



vl^it)x\it)x^^it) and 

+vl^it){x\it)x^^M + ^^itXxit) + ^AtXM 
What are the initial conditions? Again from (16) we see that 



4,(0) 



X 



(0)=0 if A,//,z/7^0 



and therefore by variation of constants we get 

x\it) xlit) 
xlit) xlit) 





-xlir) xlir) 




(63) 



(64) 



(65) 



(66) 



(67) 
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and 



x: 



X. 



■it) 
■it) 



xl(t) xl(t) 
xi{t) xi{t) 




xHt) -xI{t) 
-4(r) 



xI{t 




giir) 



(68) 



for all indices * = A/i and * = Xfiu with \, fi,iy ^ 0. But do we know the functions 
gl:{t) in all these cases? Let us look at (65). The v^{t) are all known, see (48) and 
and (51). The x\{t) have been determined in section 8, so we know all gx^it) and 
hence the a;^^(t) for A,/x 7^ from (67) and (68). These in turn give us, now for 
0, the gl^^{t) by (65) and the x^^^(t) from (67) and (68). 
It remains to determine the x\^{t) and x\^^{t) in those cases where one or several 
of the indices are zero. The values of these functions at t = T^Eq) might also be 
needed in (25) for the calculation of the partial derivatives of the Poincare map. 
The index denotes the time derivative. Knowing the Xo(f) from (59) we derive 
XQQ{t) and a;[;oo(t) as 



^00 (0 
^00 (^) 
^oo(*) 
^00 W 



X, 



000 



000 



000 



X, 



000 





m)--V2{t) by (8) 




-V22{t)y{t) 


-V222{t)y{tf + y2{t)V22{t). 



(69) 



For A,;u 7^ the x^^it) are obtained by (53), (54) and (58), (60), (61) from the 
known x\{t) and similarly the x^^^ii^) by (63), (64) and the first equation of (65) 
from the x\{t) and x\^{t). Finally, differentiating (53), (54), (58), (60) and (61) we 
see that for A G {1,2,3,4} 



x 



OOA 



X, 



OOA 



X, 



OOA 



X, 



OOA 



{t) = 
(t) = 
(t) = 



-Vu{t)x{{t) 

-V22{t)xl{t) 

-Vn2{t)y{tyxit) + Vn{t)xUt) 

-V222{t)y{t)xl{t) + V22{t)xl{t), 



(70) 



in particular x^Q^it) = for r + A odd and r, A e {1, 2, 3, 4}, and 

2^005^ = -V22{t)xl{t) + F^2{t)y{t) + Fu{t)V2{t) 

4o5it) = 

^05^ = -^222(t)|/(t)xi(t) + y22(t)4(t)-F22(0l/W-^24(t)V2(^). 



(71) 



In principle we now have all we need to calculate the partial derivatives of the 
Poincare map. 
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10 Summary of the procedure 

Once the numerical prercqiiisitcs of section 2 are established, we get the x\{t) as 
described in section 8 almost without further calculation, the only exceptions are 
and Xg, for which the integral (62) has to be evaluated. As explained in section 9, 
we also have the x\^{t) and x\^^{t) in those cases where at most one of the indices 
A, /X, V is different from zero. Next determine the x\jj:) for A, 7^ by calculating 
g'x^{t) from the first equation of (65) and applying (67) and (68). Then we get XQ)^^{t) 
for A, /i 7^ by (63), (64) and by the first equation of (65) without new integration. 
Also the fiow time derivatives Z^^ of (42) are now known. Finally, we now have 
the Qx^yif) for A,//, 7^ from the second equation of (65) and we can calculate 
the corresponding x\^^{t) as the integrals (67) and (68). Taking values at t = T(eo) 
of all these functions and applying (25), we obtain the 38 partial derivatives of the 
Poincare map at (0,0,£o)0) we wanted. 

For the computer, these instructions may be good enough, but a person might 
want to see step by step what is going on. For this we have some choice in which 
order to proceed. We will first describe all those steps that are not connected with 
the deformation question. 



11 The undeformed system 

Step 1. Choose the potential V{x, y) to be studied, with |^(0, y) = 0, choose one of 
its libration families on the y-axis and a reference point Eq for the energy parameter 
e = E — Eq. Choose a fixed Eq at which the bifurcation behavior of the libration 
shall be predicted. 

Step 2. Set up a first part of the numerical prerequisites, namely y{t), ^i(t), ^2(^) 
and their first derivatives, as described in section 2, including the period r(£o)- 
Define the four functions x\{t) on [0,T(£o)] with r, A e {1,3} by (56). 

Step 3. Collect the Jacobian matrix of the Poincare map at (0, 0, eo), or monodromy 
matrix of our librating orbit, as 

Qq Qv\^( ^im^o)) x\{T{s^)) \ 
P, P, ) \ x\{T{e,)) xl{T{eo)) ^ ^ 

according to (25), (26) and (55). If the trace Qq + Pp is different from +2, the 
fixed point is regular and the orbit will not bifurcate. In this case the procedure 
may stop here, since then we might not be interested in the higher derivatives. But 
Qq-\- Pp — 2 is not a technical necessity for going on. 

Qqe + Ppe of the tracc. 



(73) 



Step 4. Now we will calculate the £-derivative Try^(£o) 
From (25) we find 

V2 V 2 ^2 

p - -J^rV^ — 
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at (0, 0, £o)- Here we need the remaining functions and ?72(^) of the prerequisites: 
they define the x\{t) on [0, T{eo)] with r, A G {2, 4} by (57). Apart from the factor 
X2{T{eo)) in the first summand, they are needed as functions on [0, T{so)] to calculate 
and X32 in the second summand by integration (67), because the inhomogeneities 
glij,{t) and gxij,{t) for A e {1,3} and e {2,4} turn out by (65) to be 

glit) = -Vn2{t)x\{t)xl{t). 

If step 3 has shown Eq to be singular (Tr^(£o) = 2), then after completion of step 4 
we know if it is a cross-bifurcation, that is if Tr^(£o) 7^ 0. 

Step 5. Is this cross-bifurcation transcritical? To answer this question, we need 
the monodromy matrix (72) from step 3 and the second partial derivatives of Q 
and P by the variable q and p at (0, 0, Sq), that is the (x* o Z)^^ for i G {1,3} and 
m G {1,2}, to see if Pgg ^ 0, where the 'tilde' denotes adapted coordinates. By 
(26) and from our knowledge of the x\, the second equation of (25) reads 

Qqq ^11 
Qqp = x\^ 

Qpp — Xgg 

P - x' 
P — 

Ppp — 

at t = T{eo). The x\^{t) for i. A,// G {1,3} have to be calculated from (67) by 
integration with inhomogeneities 

= -Vui{t)xi{t)xl{t). 



Step 6. If the cross-bifurcation is not transcritical, then we are interested in 

o" {r\\ — ^Qm^QP ~ QpPqqq /^^x 

since the bifurcation is fork-like if and only if £^(0) ^ 0, and its sign and absolute 
value describe geometric properties of the fork. We need information beyond the 
first five steps only for P^^ and Pqqq- In the present step 6 we will take care of Pq^. 
For this, we only have to complete step 4 by the calculation of Qp^ and Pg^, which 
turn out to be ^ ^ 

^P'^ ~ TAT/ ^3-^2 "I" t7~"^32 
V'2 V2 V2 

(78) 



~ T/2V2 V^2 
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similar to (73), with 0:32 and determined by integration (67) with inhomogeneities 
given by (74). 

Step 7. To calculate Pqqq, we will need the third order partial derivatives of Q and P 
by q and p. As a preparation we will now extend steps 4 and 6 to the determination 
of all x^^(t) for r, A,/i e {1,2,3,4} by (65) and (67), (68). In all these cases we 
have 

glit) = gUt) = (79) 

and 

-Vuiit)x\{t)xl{t) if A,/iG{l,3} 
gW = { -Vu2{t)xl{t)xl{t) if A G {1,3} and //e {2,4} 

if A,Aie{2,4} 

(80) 

-Vn2{t)xl(t)xl(t) if X, lie {1,3} 
gUt) = { if Ae{l,3}and^e{2,4} 

-V222{t)xl{t)xl{t) if A,/ie{2,4} 

Note that therefore 

^Ut)-^Ut)-0 if A,Aie{2,4} and 

(81 ) 

^Ut)=4,{t)=0 if A €{1,3} and ^g{2,4}. 

To calculate the other functions a;^^(t) for r. A, e {1, 2, 3, 4} by (67) and (68) would 
be step 7. 

Step 8. It is now a suitable moment to improve upon (42) and write out the second 
order flow time derivatives as 

Vn 4 



70 _ 
^11 — 


1 4 
"^^^ 


70 _ 


1 4 
1-2"^^=^ 


^13 — 
^23 — 


1 4 

1/2"^^^ 



0, 



-x 



4 



(82) 



again at (0, 0, £05 0), with the x'^ to be taken at t = T{eo) in the notation introduced 
with (47). We will need these numbers when we apply the third equation of (25) to 
calculate Qqqq, Qqqp, • • • , Pppp- But bcfore we can do this, we have to determine the 
x\^,{T{eo)) for r. A, yu, z/ e {1, 2, 3, 4}. 



16 



Step 9. This step is parallel to and based on step 7. We read the inhomogeneities 
g\^y{t) from the second equation of (65) and then refer to (67), (68) to determine 
the x\ {t) by integration. Again 



for all A, e {1, 2, 3, 4}, and the g\^i,{t) are given in the four cases a)-d) as follows, 
a) If \,iJL,v& {1, 3}, then 



9l,M = -Fiiii(t)xi(t),Ti(t).Ti(t) 

-V,ii{t){xi{t)x'^,{t) + xl{t)xl^{f) + xl{t)xi^{t)) (84) 
-V,^2{t){x\{t)xl,{t) + xl{t)xl^{t) + xl{t)xl^{t)) , 

b) if A, /X e {1, 3} and u e {2, 4}, then 

9L.(t) = -Vrn2(t)xi{t)xl,it)xl{t) 



(83) 



Vni{t){xi{t)xl^{t) + xl{t)x^ 
Vr,2{t)xl{t)xUt) , 



it)) 



(85) 



c) if A e {1, 3} and /U, G {2, 4}, then 

9l.it) = -Vn22{t)xi{t)xl{t)xl{t) 

-Vn2(t){x\(t)xl,(t) + xl(t)xl^(t) + xl(t)x 



(86) 



d) if \ii,v e {2,4}, then 

9Lu{t) = t 



0. 



(87) 




(88) 



b) if A, e {1, 3} and v e {2, 4}, then 

9i,M = -V,,22{t)x\{t)xl{t)xl{t) 

-Vrr2{t){x\{t)xl,{t) + xl{t)x.^ 

-V222{t)xl{t)xUt) , 



(89) 



c) if A e {1, 3} and /i, u e {2, 4}, then 
9LAt) = 0, 



(90) 
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d) if X,iJ,,ue {2,4}, then 

-V222{t){xlit)x%it) + Xlit)xl,it) + Xlit)x%it)). 

Note that in particular 

= W =0 if A, ^, 1/ e {2, 4} 

= W = if A G {1, 3} and ^, G {2, 4}. 

The other functions x\^^,{t) for r. A, G {1, 2, 3, 4} have to be calculated by inte- 
grations (67) and (68) in step 9. 

Step 10. We can now write down the eight missing numbers Qqqq, Qqqp, • • • , Pppp, 
thereby completing the first part of our program, the part that is not concerned 
with deformation. The Z"^^ and with i,m G {1,2} in the formulas are taken 
from (36) and (82). The third equation of (25) gives: 

Qqqq ~ -^111 ~^ '^''^12'^11 ~^ 3x'^Z^^ 

Vg<3P — -^-113 -^-32^11 T" •^3^11 I" ^-^-1^12 

f) — _|_ ^1 72 I ™3 7O I 9^3 7O 

Vgpp — '^'133 ^ -^12^22 ^ •''1^22 ^ ^•*'1'^12 

Qppp ~ -^333 3^32^22 "I" "^^3^22 



Pqqq ~ -^111 + 3.T'^2-^11 3Viix|[^2]^ 

-Pggp = 3;5i3 + x\2Z'li — VnxlZ^^ — 2Viix\Z^2 

Pqpp — •'^133 ^" •^12'^22 ~ ^11^3'^22 ~ '^Viix\Z^2 

Pppp = X333 + 3a:|2-^22 ~ 31^11X3^22 



(93) 



at t = T(£c 



12 Deformation 

Step 11. Choose a deformation term F{x,y,Px,Py) satisfying the libration pre- 
serving condition (6). No new 'numerical prerequisites' are required, we can start 
right away calculating xl{t) and xl{t) from (62), which is the same as (68) with 
gf(t) = F^{t) and gi{t) = -F2{t). For the other xi{t) see (50) and (61). 

Step 12. Now we can determine the sixteen functions (t) for r,A G {1,2,3,4}. 
Again we derive the corresponding Qxsit) from (65). It turns out that gxsit) = if 
r -I- A is odd and hence also 

xl^(t)^0 if r + A is odd. (94) 

For A G {1,3} we get 

gl,{t) = F,,{t)xi{t) + F,,{t)xl{t) 

ghit) - -Fii(t)xi(i)-Fi3(t)xi(t)-lii2(t)xi(t)xi(t), 
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and if A e {2,4}, then 

ghit) = F24{t)xl{t) + Fu{t)xl{t) 



Step 13. Now we collect the needed flow time and starting point derivatives at 

(0, 0, So, 0) that involve the deformation parameter = S, namely the and Z^^ 
for A = 0, 2 and £ = 1, 2. From (36), (41) and (42) we get 



Z? 



F 
% 
1 



4 . 



(97) 



and 



72 _ 72 _ 

^14 — ^24 



7 _ 

^14 ~~ T/ 15 



V2 (98) 



Zo, 



24 — Tr ^35 



^2 



2 



at (0,0, £0,0). 

Step 14. We can now calculate those second derivatives of the Poincare map in 
which the deformation parameter is involved: 

Qq5 — •'^'15 + •^12-^4 ~l~ x\Z^ 

QpS — 3:35 + 3:32^4 + x^Z^ 

+ x\^Zl - V^^x\Zl 



Pq8 — X]r^ 
PpS — 3^35 + ^^2^1 — VllxlZ^ 



at t = T{eo). But before we reach the third derivatives of Q and P involving 5, we 
have to take one more step. 

Step 15. We have to determine the x\^^{T{eo)) for A,^ e {1,2,3,4}, again by 
integration (67) and (68), with inhomogeneities g\^j^^{t) as follows. 

la) If A, e {1, 3}, then 



+Fn3(t) [xi(t)xl(t) + xl(t)xl(t)] 

+F-ss3it)xlit)xlit) 
+Fi3{t)xl^{t)+F,3{t)xl^{t), 



(100) 



lb) if A e {1, 3} and e {2, 4}, then 

dUi^) = F2ss(t)xl(t)xl{t) + Fsu{t)xl(t)xl(t) 
+Fi23{t)x\{t)xl{t) + F,u{t)x\{t)xl{t) 
+Fi3it)xUt) + Fssit)xl^{t), 



Ic) if A, /X e {2,4}, then 



2a) If A,/x, e {1,3}, then 



Fn,{t)x\{t)xlit) 

+ Fi34(t) [xiit)xlit) + Xl{t)xl{t) 

+Fss4it)xUt)xlit) 
+F24{t)xl^{t) + Fu{t)xi^{t), 



2b) if A e {1, 3} and n e {2, 4}, then 



2c) if A,/x, e {2,4}, then 

gUit) 



F22,{t)xlit)xl{t) 

+F2u{t)[xl{t)xlit) + xl{t)xl{t) 

+Fui{t)xlit)xf^it) 
+F2A{t)xUt) + F,,{t)xUt). 



3a) If A, //, G (1, 3}, then 

gl,S) = -yiii2{t)x\{t)xi{t)xi{t) 

-T^iiW [a;i(i)45(i)+x;,(i)a;i5(t) 
-Ki2(t) [x\{t)x%{t) + a:;,(t)x25(t) + a:2(t)a:i^(t) 
-F,,^{t)x\{t)xl{t) 
-Fiisit) [xiit)xl{t) + xl{t)xl{t) 

-F,,,{t)xl{t)xl{t) 
-Fu{t)x\^{t)-Fn{t)xl{t), 



3b) if A e {1, 3} and ^ e {2, 4}, then 

gUit) = -V^u2(t)xl(t)xl{t)xl(t) 

-Vi22{t) [x\it)x%it) + Xlit)x\,it) + Xl{t)x\^{t) 

-Fn2{t)xi{t)xl{t) - Fu4{t)x\{t)xf,{t) 
+F^2z{t)xl{t)xl{t) + Fi34(t)a:i(t)4(t) 
-Fii{t)x\^{t)-F,^{t)xl^{t), 



3c) if A,/x e {2,4}, then 



4a) If A,^ e {1,3} then 

gUit) = -V,,22{t)xl{t)xl{t)xl{t) 

-V,,2{t)[x\{t)xlS) + ^l{t)x\S)] 

-V222{t)xl{t)xl^{t) 

-Fn2it)x\it)xlit) - F2s-iit)xlit)xl{t) 
-Fusit) [x\{t)xl{t)+xl{t)xl{t)_ 

-F22(t)xl^(t)-F2,(t)xi^(t), 

4b) if A e {1, 3} and e {2, 4}, then 

gU(t) = -y222{t)xl{t)xl,{t), 

4c) if A,// e {2,4}, then 

gUit) = -V2222{t)xl{t)xl{t)xl{t) 

-V222{t) [xl{t)x%{t) + Xl{t)xl,{t) + Xl{t)xl^{t) 
-F222{t)xl{t)xl{t) 

-F22A{t)[xl{t)xl{t) x\{t)xl{t) 
-F2u{t)x\{t)xl{t) 

-F22{t)xl^{t)-F2,{t)xi^{t). 



Note that in particular 

^UW = ^V5W = if A,//e{2,4}. 
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Step 16. This then will be the last step in the procedure to determine the 38 
numbers we set out to calculate. Only six of them remain, and these are given by 



QqqS 


~ ^115 


+ 


■^112^4 


+ 


^,3 7O 


+ 


0^3 70 


QqpS 


— ^135 


+ 


^1 72 
■^123^4 


+ 


3 7O 
-^13^4 


+ 


3 7O I 3 7O 
■^1^14 "T -('3^24 


QppS 


— ^335 


+ 


^1 72 
-^233^4 


+ 


^3 7O 
•^33^4 


+ 


9^3 7O 


PqqS 


— ^3 
~ -^115 


+ 


^3 72 
■^112^4 










PqpS 


— ^3 
~ -^135 


+ 


^3 72 
-^123^4 










PppS 


— ^3 
~ -^335 


+ 


^3 72 
■^233^4 




Viix\^Z^ 




2VnxlZl, 



(113) 



at t = r(£o); as before, or explicitly: the left hand sides are meant to be taken at 
{q,p,e,6) = (0,0,£o)0), thus on the right hand sides the x^^{t) at t = T{eo), the 
V^{x,y) at {x,y) = (0, ymax(£o))) and the are also to be taken at {q,p,e,6) = 
(0,0, £0,0), as in (41), (42), (82), (97), (98) where they have been calculated from 
xl(T{eo)), KK(0,?/max(i^o)) and various F*(0, ?/max(£o), 0, 0). 
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